package mathematics;

import org.apache.commons.math.MathException;
import org.apache.commons.math.distribution.NormalDistributionImpl;

public class FunctionNLM implements IFunction
{

	double sigmabs;
	
	public FunctionNLM(double sigmabs){
		this.sigmabs=sigmabs;
	}
	
	
	@Override
	public double getValue(double[] args) throws MathException{
		double res=0;
//		(forw(x)-1)^2 + (atmcp(x)-bscp)^2
		double x=args[0];
		double y=args[1];
		double bscp = pnorm(sigmabs * .5) - pnorm(-sigmabs * .5);
		res=Math.pow(forw(x,y)-1,2)+Math.pow((atmcp(x,y)-bscp), 2);		
		
		return res;

	}
	
	private static double forw(double x, double y) throws MathException{
		double res=0;
//			forw <- function(asw) {
//			 a <- max(asw[1],0);
//			 sw <- abs(asw[2]);
//			  wmin <- sqrt(a) - 1
//			  sw/sqrt(2*pi)*exp(-wmin^2/(2*sw^2))*(wmin + 2) +
//			  (sw^2 + 1-a)*pnorm(-wmin/sw)
//			 }
		double a = Math.max(x, 0);
		double sw = Math.abs(y);
		double wmin = Math.sqrt(a)-1;
		res=sw/Math.sqrt(2*Math.PI)*Math.exp(-Math.pow(wmin, 2)/(2*sw*sw))*(wmin+2)+(sw*sw+1-a)*pnorm(-wmin/sw);	
		return res;
	}

	private static double pnorm(double sigma) throws MathException
	{
		NormalDistributionImpl n = new NormalDistributionImpl(0.0, 1);
		return n.cumulativeProbability(sigma);
	}

	
	
	private static double atmcp(double x, double y) throws MathException{
//		 atmcp <- function(asw) {
//		 a <- asw[1]; 
//		 sw <- asw[2];
//		  wmin <- sqrt(1 + a) - 1
//		  sw/sqrt(2*pi)*exp(-wmin^2/(2*sw^2))*(wmin + 2) +
//		  (sw^2-a)*pnorm(-wmin/sw)
//		 }  
		double res=0;
		double a=x;
		double sw=y;
		double wmin=Math.sqrt(1+a)-1;
		res=sw/Math.sqrt(2*Math.PI)*Math.exp(-wmin*wmin/(2*sw*sw))*(wmin+2)+(sw*sw-a)*pnorm(-wmin/sw);
		return res;
	}
	
}
